www.gusucode.com > Heart Sound Classifier工具箱源码matlab程序代码 > Heart Sound Classifier/HeartSoundClassificationNew/HelperFunctions/dominant_frequency_features.m

    function [maxfreq, maxval, maxratio] = dominant_frequency_features(data, fs, cutoff, plotflag)
% This function estimates the dominant frequency of each channel (COLUMN) of the
% input data.  First, it calculates the power spectrum of the signal.  Then it finds the 
% maximum frequency over as well as the "energy" at that point; options below allow normalization
% by the total or mean energy.  Note this does not restrict the range of frequencies to search 
% for a maximum.
% 
% data         matrix containing data where each COLUMN corresponds to a channel
% fs           sampling rate (Hz)
% cutoff       cutoff frequency (Hz)
% maxfreq      frequency at which the max of the spectrum occurs (Hz) (ROW VECTOR)
% maxratio     ratio of the energy of the maximum to the total energy (ROW VECTOR)
%Copyright (c) 2016, MathWorks, Inc. 

[~, nchannels] = size(data);

% nfft = 2^nextpow2(npoints*padfactor);
nfft = 4096;
f = fs/2*linspace(0,1,nfft/2);
cutoffidx = length(find(f <= cutoff));
f = f(1:cutoffidx);

% calculate the power spectrum using FFT method
datafft = fft(data,nfft);
ps = real(datafft.*conj(datafft));

% keep only the non-redundant portion
ps = ps(1:cutoffidx,:);
for ich = 1:nchannels
    ps(:,ich) = ps(:,ich)/sum(ps(:,ich));
end

if(plotflag)
    figure, 
    subplot(211), plot(data)
    subplot(212), plot(f,ps);
    axis([0 cutoff 0 max(ps)]);
    pause
end

% locate max value below cutoff
[maxval, maxind] = max(ps);
maxfreq = f(maxind);

% calculate peak energy by summing energy from maxfreq-delta to maxfreq+delta
% then normalize by total energy below cutoff
delta = 5;  % Hz
maxratio = zeros(1,nchannels);
for ich = 1:nchannels
    maxrange = f>=maxfreq(ich)-delta & f<=maxfreq(ich)+delta;
    maxratio(ich) = sum(ps(maxrange,ich)) / sum(ps(:,ich));
end